%% Uncertainty of Parameter Estimates I

%We will use a 'Monte-Carlo' method to estimate the error in our parameter
%estimate

%Setup model:

%1st order polynomial
a=1;
b=5;

x=(0:100);
y=a+b*x;

%Gaussian (normal) noise
sigma=20;

%Generate noisy data
y_N=normrnd(y,sigma);

%Setup fminsearch and find optimum values. 
Options=optimset();
Theta0=[min(y_N); (y_N(end)-y_N(1))/(x(end)-x(1))];
[Theta_Found]=fminsearch(@firstOPoly_Norm,Theta0,Options,x,y_N);

Afound=Theta_Found(1);
Bfound=Theta_Found(2);

%If Afound and Bfound are correct, what would the distribution of 
%parameters estimates look like with different noise realizations?

%Do 'Monte-Carlo' simulations to find distribution of 'a' and 'b'
y=Afound+Bfound*x;

Ntrials=1000;

%Setup vectors to store result of the trials
Avec=ones(Ntrials,1);
Bvec=ones(Ntrials,1);

for nn=1:Ntrials
    
    %Generate a set of noisy data
    y_N=normrnd(y,sigma);
    
    %Estimate 'a' and 'b' from this data set
    Theta0=[min(y_N); (y_N(end)-y_N(1))/(x(end)-x(1))];
    [Theta_Found]=fminsearch(@firstOPoly_Norm,Theta0,Options,x,y_N);
    
    %Store in vector
    Avec(nn)=Theta_Found(1);
    Bvec(nn)=Theta_Found(2);

end

%Look at result with histograms
figure
hist(Avec,20)
xlabel('A')
ylabel('Counts')

figure
hist(Bvec,20)
xlabel('B')
ylabel('Counts')

%Notice they are Gaussian shaped.  We can then give the 'standard error'
%as the standard deviation of the parameter estimates.

fprintf('We found: a = %g +/- %g\n',Afound,std(Avec))
fprintf('We found: b = %g +/- %g\n',Bfound,std(Bvec))







